%%%%%%%%%%%%%%%%%%%%%%%%%%% asme2ej.tex %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Template for producing ASME-format journal articles using LaTeX    %
% Written by   Harry H. Cheng, Professor and Director                %
%              Integration Engineering Laboratory                    %
%              Department of Mechanical and Aeronautical Engineering %
%              University of California                              %
%              Davis, CA 95616                                       %
%              Tel: (530) 752-5020 (office)                          %
%                   (530) 752-1028 (lab)                             %
%              Fax: (530) 752-4158                                   %
%              Email: hhcheng@ucdavis.edu                            %
%              WWW:   http://iel.ucdavis.edu/people/cheng.html       %
%              May 7, 1994                                           %
% Modified: February 16, 2001 by Harry H. Cheng                      %
% Modified: January  01, 2003 by Geoffrey R. Shiflett                %
% Use at your own risk, send complaints to /dev/null                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% use twocolumn and 10pt options with the asme2ej format
\documentclass[twocolumn,10pt]{asme2ej}

%\usepackage{psfig} %% for loading postscript figures

% WILLIAM: acrescentei packages abaixo
\usepackage{tikz}
\usepackage[tight]{subfigure}
\usepackage[ugly]{units}

%% The class has several options
%  onecolumn/twocolumn - format for one or two columns per page
%  10pt/11pt/12pt - use 10, 11, or 12 point font
%  oneside/twoside - format for oneside/twosided printing
%  final/draft - format for final/draft copy
%  cleanfoot - take out copyright info in footer leave page number
%  cleanhead - take out the conference banner on the title page
%  titlepage/notitlepage - put in titlepage or leave out titlepage
%  
%% The default is oneside, onecolumn, 10pt, final


\title{A 3D Adaptive Mesh Generation Approach using Geometric Modeling with Multi-Regions and Parametric Surfaces}

%%% first author
\author{Antonio C. O. Miranda
    \affiliation{
	Professor\\
	Department of Civil Engineering\\
	University of Bras\'ilia\\
	Bras\'ilia - Brazil\\
    Email: acmiranda@unb.br
    }	
}

%%% second author
%%\author{Rafael A. Sousa
%%    \affiliation{
%%	Graduate Student\\
%%	Department of Civil Engineering\\
%%	Pontifical Catholic University of Rio de Janeiro\\
%%	Rio de Janeiro - Brazil\\
%%    Email: rfldesousa@hotmail.com
%%    }	
%%}

%%% third author
\author{William W. M. Lira\thanks{Campus A. C. Sim\~oes - Av. Lourival Melo Mota, Cidade Universit\'aria - Macei\'o - AL, Brazil. CEP:57072-900. Tel.: +55-82-3214-1304; fax: +55-82-3214-1303}
    \affiliation{
	Professor\\
	LCCV - Laboratory of Scientific Computing\\
	and Visualization\\
	Technology Center\\
	Federal University of Alagoas\\
	Macei\'o - Brazil\\
    Email: william@lccv.ufal.br
    }	
}

%%% fourth author
\author{Joaquim B. Cavalcante-Neto
    \affiliation{
%%    \\
	Professor\\
	Department of Computing\\
	Federal University of Cear\'a\\
	Fortaleza - Brazil\\
    Email: joaquim@lia.ufc.br
    }	
}

%%% fifth author
\author{Luiz F. Martha
    \affiliation{
	Professor\\
	Tecgraf - Computer Graphics Technology Group\\
	Department of Civil Engineering\\
	Pontifical Catholic University of Rio de Janeiro\\
	Rio de Janeiro - Brazil\\
    Email: lfm@tecgraf.puc-rio.br
    }	
}

\begin{document}

\maketitle    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
{\it This work presents a methodology for adaptive generation of 3D finite element meshes using geometric modeling with multi-regions and parametric surfaces, considering a geometric model described by curves, surfaces, and volumes. This methodology is applied in the simulation of stress analysis of solid structures using a displacement-based finite element method, and may be extended to other types of 3D finite element simulation. The adaptive strategy is based on an independent and hierarchical refinement of curves, surfaces, and volumes. From an initial model, new sizes of elements obtained from a discretization error analysis and from geometric restrictions are stored in a global background structure, a recursive spatial composition represented by an octree. Based on this background structure, the model's curves are initially refined using a binary partition algorithm. Curve discretization is then used as input for the refinement of adjacent surfaces. Surface discretization also employs the background octree-based refinement, which is coupled to an advancing front technique in the surface's parametric space to generate an unstructured triangulated mesh. Surface meshes are finally used as input for the refinement of adjacent volumetric domains, which also uses an advancing front technique, but in 3D space. In all stages of the adaptive strategy, the refinement of curves, surface meshes, and solid meshes is based on estimated discretization errors associated with the mesh of the previous step in the adaptive process. In addition, curve and surface refinement takes curve information into account. Numerical examples of simulation of engineering problems are presented in order to validate the methodology proposed in this work.
}
\end{abstract}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{nomenclature}

\entry{$\xi_i$}{discretization error ratio.}
\entry{$\xi_{avg}$}{average error ratio of finite element model.}
\entry{$\xi_{dev}$}{standard deviation error ratio of finite element model.}
\entry{$\nu$}{Poisson's coefficient of material.}
\entry{$\theta_{ij}$}{estimated angle between curve tangents between vertices $i$ and $j$.}
\entry{$\theta_{max}$}{maximum allowed curvature angle of curves.}
\entry{$\eta$}{relative percentage error in the energy norm.}
\entry{$\eta_es$}{estimated relative percentage error.}
\entry{$\eta_max$}{maximum allowable relative percentage error.}
\entry{$\vec{\sigma_{U}}(U, V)$}{gradient vector in direction $U$ at a parametric point $(U,V)$ on a surface.}
\entry{$\vec{\sigma_{V}}(U, V)$}{gradient vector in direction $V$ at a parametric point $(U,V)$ on a surface.}
\entry{$\sigma$}{exact stress fields.}
\entry{$\hat{\sigma}$}{finite element stress fields.}
\entry{$\bar{\sigma}$}{recovery stress fields.}
\entry{$\vec{\sigma_{V}}(U, V)$}{gradient vectors, direction $V$, at a parametric point $(U,V)$ on the surface.}
\entry{$\Omega$} {the problem definition domain.}
\entry{$\vec{N}_{2D}$}{unit vector between points $m_{AB}$ and $N_{C}$.}
\entry{$ch_i$}{characteristic size of a finite element.}
\entry{$d_{ij}$}{distance between vertices $i$ and $j$ on a curve.}
\entry{$D$}{constitutive matrix.}
\entry{$||e||$}{energy norm for the displacement error.}
\entry{$||e||_{es}$}{estimated energy norm for the displacement error.}
\entry{$E$}{Young modulus of material.}
\entry{$F, G, K$}{metrics of a point on surface.}
\entry{$|J|$}{determinant of the Jacobian transformation matrix.}
\entry{$h_{2D}$}{height of triangle in parametric space.}
\entry{$h_{3D}$}{height of triangle in three-dimensional space.}
\entry{$H_{new,old}$}{new and old quadtree cell sizes.}
\entry{$L_i$}{length of an edge.}
\entry{$m_{AB}$}{midpoint of edge $AB$.}
\entry{$N_{A,B}$}{normal vectors at points $A$ and $B$ on surface.}
\entry{$N_{C}$}{vertex to create a triangle with edge $AB$.}
\entry{$p_{i,j}$}{vertices $i$ and $j$ on a curve.}
\entry{$p_{x,y,z}$}{pressure in  directions x, y, z.}
\entry{$s_{ij}$}{circle arc from vertices $i$ to $j$.}
\entry{$T_{i,j}$}{tangent vectors at vertices $i$ and $j$ on a curve.}
\entry{$U, V$}{parametric coordinate system.}
\entry{$||u||$}{square root of twice the strain energy.}
\entry{$||\hat{u}||$}{energy norm obtained from the finite element solution.}
\entry{$u, v, w$}{displacement in directions $x, y, z$.}
\entry{$Vol_i$}{volume of a finite element.}
\entry{$x, y, z$}{Cartesian coordinate system.}
\end{nomenclature}

\noindent Keywords: 3D Adaptive Analysis, Finite Element Method, Mesh Generation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}


In numerical simulations using the Finite Element Method (FEM), two important aspects to be considered are the automatic generation of the model's finite element mesh and the definition of the level of refinement associated with this mesh. For the first aspect, there are a variety of algorithms with different techniques to generate planar, surface, and volumetric meshes. The second aspect, i.e., the level of refinement, is usually defined manually by a specialist based on his/her own experience. However, this refinement should consider the fact that the density of the generated elements varies according to the gradient of the obtained solution, which is initially unknown. In this context, this work presents a methodology to automate the refinement process of three-dimensional meshes based on adaptive technique procedures. In this work, the proposed methodology is applied to stress analysis of solid structures using a displacement-based finite element formulation. However, since the methodology essentially treats geometric modeling and mesh generation aspects of the problem, it could be used in other types of 3D finite element simulation.

Adaptive procedures try to automatically refine and coarsen a mesh, relocate its nodes, or adjust its cells to improve response accuracy. Usually, the computation begins with a trial solution obtained from a coarse mesh. The discretization error of this solution is estimated. If it fails to satisfy a prescribed accuracy metric, adjustments are made to achieve the desired solution with minimal effort. Common procedures are \cite{Mackerle2001, McRae2000}: local/global refinement and/or coarsening of a mesh (\textit{h-refinement}), relocating or moving a mesh (\textit{r-refinement}), and locally varying the polynomial degree of elements (\textit{p-refinement}).  

Some strategies have been proposed to efficiently automate the 3D mesh refinement process. These strategies can be divided in two approaches: local and global refinement. In local refinement, the process uses an initial mesh and locally, using a set of elements, refines or coarsens elements in the mesh. Most works in the literature are based on this approach, as described ahead. Kallinderis and Vijayant \cite{Kallinderis1993} and Muthukrishnan \cite{Muthukrishnan1995} present an adaptive grid scheme based on the division/deletion of tetrahedral cells. Golias and Tsiboukis \cite{Golias1994} and Golias and Dutton \cite{Golias1997} employ a set of topological Delaunay transformations of tetrahedral elements and a technique for node reposition. Lee and Lo \cite{Le1997,Le1999} approach mesh refinement by inserting additional nodes at the midpoint of the longest or quasi-longest line segment of the mesh that bisects the original edges to generate new elements. In Merrouche \cite{Merrouche1998}, the mesh adaptation is achieved by a 3D bisection method. De Cougny and Shephard \cite{Cougny1999} present an adaptive scheme based on subdivision patterns (for refinement), edge collapsing (for coarsening), and mesh optimization (following refinement and coarsening). Lee et al. [9] only increase the order and density of 3D finite element meshes. Lee and Xu \cite{Lee2005} generate a surface mesh for the mid-surface of the thin-walled structure, controlling element size, and convert the surface mesh to a 3D solid mesh by extrusion. More recently, Zhang et al. \cite{QianZhang2012, Zhang2010405} generates tetrahedral and hexahedral meshing in multi-material domains using grid-based method that employs an octree structure, refining meshes also locally. 

In global refinement, on the other hand, at each refinement step, the entire mesh is deleted and another is generated based on new sizes of elements obtained from a discretization error estimation analysis. This process is used by Kettil et al. \cite{Kettil2003} only in regions with complicated parts. Hughes et al. \cite{Hughes2005} refine the structured meshes with NURBS surfaces. Our work presents a methodology that employs the global h-refinement approach.

A previous paper \cite{Cavalcante-Neto1998} proposed a two-dimensional self-adaptive strategy that was able to perform simulations involving automatic generation of meshes and adaptive methods. Other works have considered the same problem through different approaches, such as the study by Mark Shephard's team \cite{Baehmann1989}, from the Rensselaer Polytechnic Institute. Cavalcante-Neto \cite{Cavalcante-Neto1998b} proposed a technique for the generation of volumetric meshes of tetrahedral elements for arbitrary region domains. Combining this technique with the implementation of 3D error estimators, the authors defined a prototype of an environment for adaptive generation in three dimensions. However, in this previous work, the complete process was not performed automatically, i.e., the mesh was generated independently from the error estimation and had to be manually combined. Moreover, despite treating multi-regions, it was not very efficient, and parametric surfaces, which are used in several types of simulations, were not considered.

This paper aims to present a methodology for adaptive generation of three-dimensional finite element meshes, using geometric modeling with multi-regions and parametric surfaces. Basically, the whole mesh adaptive process involves three steps: (1) analysis of a finite element (FE) model with discretization error estimation; (2) construction of a background structure to store new FE sizes that take into account te estimated discretization error and curve and surface curvature; and (3) hierarchical refinement of a FE model that is represented geometrically by curves, surfaces, and volume regions. This process may be repeated until a desired maximum allowed error metric is achieved. The methodology described herein covers in detail only the last two steps, since discretization error estimation can be computed through different processes \cite{Mackerle2001}. However, for the sake of completeness, a summary of the adopted discretization error estimation technique is presented.

The paper is organized as follows. In the next section, some background on the adopted discretization error estimation technique is summarized.  Section~\ref{sec:RefinementStrategy} explains in detail the proposed adaptive refinement strategy. The following section describes all the steps required to generate the background data structure that is used to define FE sizes in the adaptive process. Section~\ref{sec:HierarchicalRefinement} presents the hierarchical refinement of curves, surfaces, and volumes. Five examples of adaptive refinement are presented in Section~\ref{sec:Examples} to evaluate result convergence and performance and the quality of the generated elements. Finally, in Section~\ref{sec:Conclusion}, the authors present comments and conclusions of this work.


\section{Theoretical Background: Discretization Error Estimation and Adaptive Refinement}\label{sec:ErrorEstimation}

This theoretical background is based on a standard technique used in the literature \cite{Zienkiewicz2006}. There are different approaches for this technique, as presented by Mackerle \cite{Mackerle2001}, and one of them is described here. This section is a summary of an existing technique and the authors do not aim to discuss it in detail. 

To compute the local error correctly, mathematical norms are introduced to measure the discretization error. The discretization error in a displacement-based finite element solution is often quantified on the basis of a norm for the strain energy error, $||e||$, which can be expressed, in terms of stresses, as

\begin{equation}
||e||^2 = \int_{\Omega} (\sigma-\hat{\sigma})^T D^{-1} (\sigma-\hat{\sigma}) d \Omega
\label{EQ_error},
\end{equation}

\noindent where $\sigma$ and $\hat{\sigma}$ are the exact and the finite element stress fields, $D$ is the constitutive matrix, and $\Omega$ is the problem definition domain.

The basic purpose of error estimators is to substitute the field $\bar{\sigma}$, obtained by means of recovery procedures (e.g. ZZ, SPR, or REP) \cite{Zienkiewicz2006, Zienkiewicz1992a, Zienkiewicz1992b} for the field $\sigma$, which is generally unknown. Therefore, the expression for computing the approximate (estimated) relative error distribution $||e||_{es}$ can be expressed as

\begin{equation}
||e||^2_{es}= \int_{\Omega} (\bar{\sigma}-\hat{\sigma})^T D^{-1} (\bar{\sigma}-\hat{\sigma}) d \Omega
\label{EQ_error_est}.
\end{equation}

Taking the finite element discretization into account and considering a specific finite element $i$, Equation \ref{EQ_error_est} can be rewritten as

\begin{equation}
(||e||^{i}_{es})^2= \int_{\Omega_i} (\bar{\sigma}-\hat{\sigma})^T D^{-1} (\bar{\sigma}-\hat{\sigma}) |J| d \Omega_i
\label{EQ_error_est_elem}.
\end{equation}

\noindent where standard isoparametric elements are assumed, $|J|$ is the determinant of the Jacobian transformation matrix, and $\Omega_i$ is the element domain. 

The energy norm for the error can be evaluated over the whole domain or part of it. The contribution of all the elements $m$ in the mesh is given by 

\begin{equation}
||e||^{2} = \sum_{i=1}^{m} (||e||^{i})^2
\label{EQ_energy_norm}.
\end{equation}

The relative percentage error in the energy norm, $\eta$, for the whole domain or part of it, can be obtained as $ \eta = ||e|| / ||u||$, where $||u||$ is the square root of twice the strain energy, and is given by

\begin{equation}
||u||= \left( \int_{\Omega_i} \sigma^T D^{-1} \sigma d \Omega_i \right)^{0.5}
\label{EQ_error_strain_energy}.
\end{equation}

A simple criterion to achieve a solution error with an acceptable level for the whole domain can be stated as $\eta_{es} \leq \eta_{max}$, where $\eta_{max}$ is the maximum allowable error defined by the analyst and the estimated error, $\eta_{es}$, is given by

\begin{equation}
\eta_{es} =  \frac{||e||_{es}}{ \left( ||\hat{u}||^2 + ||e||_{es}^2 \right)^{1/2}}
\label{EQ_error_estimated},
\end{equation}

\noindent where $||\hat{u}||$ is the energy norm obtained from the finite element solution.

A sound criterion for an ``optimal mesh'' consists of requiring that the energy norm error be equidistributed among the elements, because it leads to meshes with high convergence rates \cite{Zienkiewicz2006}. Thus, for each element $i$,

\begin{equation}
||e||^i_{es} < \eta_{max} \left( \frac{||\hat{u}||^2 + ||e||_{es}^2}{m}   \right)^{1/2} = \bar{e}_m
\label{EQ_error_medium}.
\end{equation}

\noindent By defining the ratio

\begin{equation}
\xi_i =  \frac{||e||^i_{es}}{\bar{e}_m} 
\label{EQ_error_ratio},
\end{equation}

\noindent then refinement is needed if $\xi_i > 1$. By assuming a certain rate of convergence, the value of the discretization error ratio $\xi_i$ can be used to decide the new size of the element. Thus

\begin{equation}
ch_i =  \frac{h_i}{\xi_i^{1/p}} 
\label{EQ_error_new_size},
\end{equation}
 
\noindent where $h_i$ is the initial size of the element, $ch_i$ is the desired  characteristic size of the element, and $p$ is the polynomial order of the finite element approximation.
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Adaptive Refinement Strategy}\label{sec:RefinementStrategy}

The three-dimensional geometric model has a topological description of the vertices, curves, surfaces, and regions, as well as an associated geometric description, which consists of the coordinates of the vertices and the mathematical representation of the curves and surfaces. The geometric model can contain many regions. In this environment, the attributes of the simulation, such as the properties of the materials, loads and restrictions, are associated with the geometric entities. In this framework, the entities of finite element mesh (nodes and elements) automatically receive the attributes of the geometric entities that are related to. Using this approach, it is possible to create new meshes without losing the attributes.

Figure \ref{FIG_proposed_process} illustrates the automatic adaptive strategy of the proposed refinement process. The input data are the initial volumetric mesh of problem in question and the geometric entities (curves and surfaces), and as well as their associated attributes. Initially, this mesh is numerically analyzed, the information required to initiate the adaptive procedure. Such information basically consists of numerical discretization errors associated with each volumetric element of the mesh. From these errors, the need for adaptive refinement is verified. If the results converge, the adaptive process is concluded with a final discretization.

%% FIGURE
\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-1.png}}
\caption{The proposed adaptive refinement process.}
\label{FIG_proposed_process}
\end{figure}

If convergence is not reached, the sizes of the new elements are computed based on the estimated discretization error. All the resized data are stored in an auxiliary background structure. Although many background structures are published in the literature, as reviewed by Quadros \cite{Quadros2010}, the present work uses a background octree structure, which has the advantage of not only allowing fast search procedures down to internal leaves but also representing the desired size of the elements defined by the size of the internal leaves. For these reasons, an octree is used to support the discretization of curves, surface meshes, and the volumetric mesh. 

In addition to discretization error estimation, curve and surface discretization is also required, especially when the curves and surfaces present high curvatures. In such locations, the meshes should be locally refined. Therefore, new element size data, based on the geometric information of the curves and surfaces, are computed and stored in the background octree. After this procedure, the background octree is internally finalized to provide a better transition between regions with elements of highly varying sizes.

Using the size information from the background octree, the next step consists of a three level hierarchical approach to create a new volumetric mesh. First, the curves are refined based on the size of the elements stored in the octree structure. This refinement subdivides the curves into segments with sizes consistent with those of the discretization error analysis and geometric criteria. After refining the curves, the meshes associated with each of the model's surfaces are discretized using an advancing front scheme in parametric space. This meshing scheme starts by subdividing curves on the boundary of each surface. Geometric curvature information is considered in surface refinement because the background octree takes this information into account. The last stage of the adaptive refinement process is related to the discretization of the domains of the model's regions. Such discretization uses a 3D advancing front technique that starts from the triangulated meshes associated with the boundary surfaces of each 3D region, also considering the sizes of the elements provided by the background octree. As can be seen, this adaptive meshing methodology supports multi-regions in a consistent manner, considering curve and surface curvature information in addition to the estimated discretization error.

Finally, a new discretization error analysis is performed to assess the quality of the results. If convergence is not obtained, the whole adaptive process is repeated as described above. 

The next sections detail the proposed refinement strategy.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The Background Octree}\label{sec:BackgroundOctree}

An octree is a tree data structure based on a cell with eight children. Each cell of an octree represents a cube in the physical space. Each child represents one octant of its parent. On the leaves of the tree are the computational cells of the grid. In this work, the background octree has two main objectives. The first is to develop local guidelines used to define the discretization of curves and surfaces. The second is to define the sizes of tetrahedral elements to be generated during the advancing front procedure. The octree generation includes four steps. In the first step, the octree is initialized based on the input mesh data, which are the new element sizes obtained in the discretization error analysis. The second step refines the octree based on the geometric curvatures of curves. In the third step, the octree is refined based on the geometric curvatures of surfaces. Finally, the octree is refined in order to obtain a better transition between the sizes of the elements generated in the advancing front surface refinement.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Background Octree Refinement Based on Error Analysis}\label{sec:ErrorBackgroundOctree}

Initially, a bounding cube is created based on the maximum range of the three Cartesian coordinates of the input model. This cube is the octree's root cell. In the first step of the octree refinement, each discretization error result of an element is used to determine the local depth of the subdivision. The characteristic size $ch_i$ of each element is calculated and the octree cell containing the element's central point is determined. If the size of the cell edge is larger than the calculated characteristic size, then this cell is subdivided into eight smaller cells. This process is repeated recursively and finishes when the size of the cell is smaller than the given size. This process is repeated for every element of the current FE mesh.

The characteristic size $ch_i$, which is used in the refinement of the background octree, is calculated as follows. First, the volume of the existing element, $Vol_i$, is obtained. Then, the length of edge $L_i$, considering an equilateral tetrahedron with the same volume, is estimated as


\begin{equation}
L_i = \sqrt[3]{\frac{12 Vol_i}{\sqrt{2}}}
\label{EQ_lenedge}.
\end{equation}

The $L_i$ value is used to obtain the characteristic size $ch_i$, Equation \ref{EQ_error_new_size}, as proposed by Zienkiewicz and Taylor \cite{Zienkiewicz2006}.

The background octree works as a density function to guide the adaptive process. It could be replaced by other functions if desired, although it has the advantage of also allowing fast search procedures down to internal leaves. It could also have a different orientation to better adapt to models that are not parallel to Cartesian coordinates. However, the bounding cube parallel to the Cartesian coordinates is easier and faster to implement and usually gives very good results.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Octree Refinement Based on Curve Curvature}\label{sec:CurvesOctree}

In some cases, when only the discretization error is considered in the adaptive process the new generated mesh (in the following step of the process) does not respect the actual geometry of the model's curves. This behavior occurs when parts of a curve, for example, present high curvatures in a region where the discretization error is low. In these situations, it is necessary to refine the background octree based on the curvatures of the curves to preserve the original geometric characteristics of the model.

The methodology used to refine the curves of the model based on their curvatures is a one-dimensional version of the procedure applied to discretize the background octree. The refinement of each curve employs a recursive spatial numbering technique similar to a binary tree data structure \cite{Knuth1997}, as illustrated in Figure \ref{FIG_curve_refinement_a}. 

The main purpose is to generate a discretization on a curve according to its curvatures. The curvatures are calculated for specific curve segments. At first, the whole length of the curve is considered as the segment to be tested. If the curvature of the segment is lower than the maximum allowed curvature, $\theta_{max}$, the process is interrupted. Otherwise, the segment is recursively subdivided in two segments, and each one is tested in the same way, until the maximum curvature criterion is satisfied. At the end of this process, all the curve segment sizes and their middle points are transported to the background octree, using the same procedure explained in Section~\ref{sec:ErrorBackgroundOctree}. Figure \ref{FIG_curve_refinement_a} demonstrates this curve refinement process. In this figure, a binary tree is used in each level of the recursive refinement.


\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-5.png}}
\caption{Curve refinement based on curve curvature.}
\label{FIG_curve_refinement_a}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=1.5in]{./figure/Figure-6.png}}
\caption{Approximating the curvature by circular arcs.}
\label{FIG_curve_refinement_b}
\end{figure}

In fact, the curvatures of the curves are measured indirectly considering local approximation by circular arcs. Let $T_i$ and $T_j$ be the tangent vectors at vertices $p_i$ and $p_j$ (see Figure \ref{FIG_curve_refinement_b}). Consider the circle arc $s_{ij}$ going from $p_i$ to $p_j$, $l_{ij}$ is the actual curve length from $p_i$ to $p_j$, and $d_{ij}$ is the distance between the two vertices. Then the angle $\theta_{ij}$ of $s_{ij}$ can be estimated as



\begin{equation}
\theta_{ij} \approx \arccos{\frac{T_i\cdot{T_j}}{\|T_i\|\|T_j\|}}
\label{EQ_curvature}.
\end{equation}

However, this curvature approximation is not sufficient to detect high curve curvatures.  For this reason, Equation \ref{EQ_curvature} is used together with another equation:

\begin{equation}
\frac{l_{ij}}{d_{ij}}\leq\frac{\theta_{max}}{2\cdot\sin({\frac{\theta_{max}}{2}})}
\label{EQ_curvature_approx}.
\end{equation}

To understand the need for Equation \ref{EQ_curvature_approx}, see for example the case of Level 0 in Figure \ref{FIG_curve_refinement_a} with Equation \ref{EQ_curvature}. In this case, the tangent vectors from only two vertices can be lower than the maximum allowed curvature. Considering also Equation \ref{EQ_curvature_approx}, the segment is subdivided into two segments, as shown in Level 1. In this example, the process is repeated until Level 5 to satisfy both equations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Octree Refinement Based on Surface Curvature}\label{sec:SurfacesOctree}

After the background octree is refined considering curve curvatures, the octree refinement is increased based on surface curvatures. This step captures high curvatures of surfaces, computes the required element sizes and their locations, and passes this information to the background octree. The reason to perform this step is the same as the previous stage: to represent the original geometry of the model.

As in Section~\ref{sec:CurvesOctree}, the methodology applied to refine the surfaces of the model based on their curvatures is a two-dimensional version of the procedure used to discretize the background octree: a background quadtree, which is created similarly to the one presented by Miranda and Martha \cite{Miranda2002}. This way of computing the curvatures has shown to be efficient and robust, and this is the main reason for its adoption. The background quadtree generation follows some steps:

\begin{itemize}
\item[$\bullet$] Quadtree initialization based on given boundary edges;
\item[$\bullet$] Refinement to force maximum cell size;
\item[$\bullet$] Refinement to provide minimum size disparity for adjacent cells;
\item[$\bullet$] Refinement to force minimum curvature difference between adjacent cells: this stage is explained ahead.
\end{itemize}
	
As described in detail by Miranda and Martha \cite{Miranda2002}, the first step has some modifications in relation to the original 3D algorithm \cite{Cavalcante-Neto2001}. The second and third steps have not changed. The fourth step was added to take high surface curvatures into account.
	
The fourth step of the quadtree generation refines this auxiliary structure to force a minimum curvature difference between adjacent cells. Initially, the algorithm stores in each cell gradient vectors of the quadtree evaluated at the center of the cell. Then, it computes a vector normal to the surface of each cell. Finally, the algorithm obtains the cosine of the angle between the normal vectors, $N_A$ and $N_B$, of the two adjacent cells:

%%
\begin{equation}
\cos (\theta) = \frac{N_A\cdot{N_B}}{\|N_A\|\|N_B\|}
\label{cos_theta_curvature},
\end{equation}


\noindent and compares it to a minimum value, $\cos(\theta_{min})$. This kind of comparison is similar to comparing the angle between the normal vectors and the maximum angle. If $\cos(\theta)$ is smaller than $\cos(\theta_{min})$, then a new cell size, $H_{new}$, is obtained from the current size, $H_{old}$, as

%%
\begin{equation}
H_{new} = \frac{H_{old}}{\cos(\theta_{min})}\cdot\cos(\theta)
\label{EQ_Hnew}.
\end{equation}


This new size is used to locally refine the adjacent cells of the quadtree. This process is repeated recursively for every cell. The new element sizes stored in the auxiliary surface quadtree are transferred to the global background octree. At the end of this step, the background octree is refined considering the geometric curvatures of all of the model's surfaces.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Octree Final Refinement}\label{sec:FinalOctree}

The previous step can leave large octree cells in the interior of a 3D region. In the first step of this final stage, the octree is refined to guarantee that no cell in the interior is bigger than the largest cell on the boundary. This will avoid excessively large elements in the domain interior. The octree is subsequently processed to force a single difference level between neighboring cells. This leads to a natural transition between regions with different degrees of refinement. This refinement is performed by traversing the octree and examining the difference in tree depth between adjacent cells. If the difference is larger than one level, the adequate cells are refined until the criterion is satisfied. 

Figures \ref{FIG_artificial_model_a} to \ref{FIG_artificial_model_e} depict the external appearance of the background octree of a hypothetical model following the steps described previously.

\begin{figure} 
\centerline{\includegraphics[width=2.5in]{./figure/Figure-3.png}}
\caption{A hypothetical model under applied $\sigma$ tension to explain the steps of octree construction.}
\label{FIG_artificial_model_a}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2in]{./figure/Figure-4.png}}
\caption{External appearance of background octree based on discretization error analysis.}
\label{FIG_artificial_model_b}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2in]{./figure/Figure-7.png}}
\caption{Refinement after considering curve curvatures.}
\label{FIG_artificial_model_c}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2in]{./figure/Figure-9.png}}
\caption{Refinement after considering surface curvatures.}
\label{FIG_artificial_model_d}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2in]{./figure/Figure-10.png}}
\caption{Refinement after considering maximum cell size at boundary cells and maximum difference of one level between adjacent cells.}
\label{FIG_artificial_model_e}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Hierarchical Refinement}\label{sec:HierarchicalRefinement}

After the construction of the background octree, considering the new element sizes based on the discretization error and on curvatures of the geometric model, the final step is to re-generate the mesh of the whole model. As mentioned previously, it is assumed that the three-dimensional geometric model has a topological description of the vertices, curves, surfaces, and regions, as well as an associated geometric description, which consists of the coordinates of the vertices and the  mathematical representation of the curves and surfaces. The geometric model can contain many regions. Mesh re-generation employs a hierarchical refinement of (1) curves, (2) surfaces, (3) and regions. The process is described in the following subsections.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Curve Refinement}\label{sec:CurveRefinement}

The methodology used to refine the model's curves is similar to the one presented in Section~\ref{sec:CurvesOctree}. The curve refinement process is explained by means of an example shown in Figure \ref{FIG_curve_refinement_tree}, which also depicts a binary tree data structure for each refinement level. 

At the beginning of the process, a curve is defined by its mathematical geometric description and by two nodes (initial and final points), as shown in Figure \ref{FIG_curve_refinement_tree}, Level 0. Then, the curve length and middle node are obtained. From the middle node, one can determine the cell in the background octree where this node is located. A comparison is then made to verify whether the segment size is smaller than that of the corresponding cell. If the criterion is satisfied, the curve refinement process ends considering the nodes generated so far. Otherwise, the new node is inserted (Level 1) on the curve, this curve is subdivided in two partitions (Level 2), and each one is tested in the same way, until the criterion is met (Levels 2 to 4).

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-11.png}}
\caption{Curve refinement and its corresponding binary tree.}
\label{FIG_curve_refinement_tree}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Surface Refinement}\label{sec:SurfaceRefinement}

Surface mesh generation is based on the algorithm presented by Miranda and Martha \cite{Miranda2002}. This algorithm is applied to the generation of triangular meshes on each surface with arbitrary geometry, using its parametric description. The parametric description is used because it is common and efficient, since the surface mesh is generated using two-dimensional triangulation techniques. However, additional length and angle corrections are needed to consider metric distortions between parametric and 3D Cartesian spaces. With this procedure, generated triangles present good shape in 3D space.

\begin{figure} 
\centerline{\includegraphics[width=2in]{./figure/measurements_in_2d.png}}
\caption{Measurements in parametric space.}
\label{measurements_in_2d}
\end{figure}

Figure \ref{measurements_in_2d} shows an edge $AB$ that will create a triangle with point $N_{C}$ in parametric space. $m_{AB}$ is the midpoint of points $A$ and $B$, $h_{2D}$ is the height of triangle, and $\vec{N}_{2D}$ is a unit vector from points $m_{AB}$ and $N_{C}$. This vector is normal to edge $AB$ in 3D space. Given a desired size $h_{3D}$ in 3D space obtained from the background structure, the distance $h_{2D}$ in parametric space is given by:


\begin{equation}
h_{2D} = \frac{h_{3D}}
          {
           \sqrt{
            F_{m_{AB}} U^{2}_{N_{2D}} + 2G_{m_{AB}} U_{N_{2D}} V_{N_{2D}} + K_{m_{AB}} V^{2}_{N_{2D}}
            }
          },
\label{EQ_h3d_to_h2d}
\end{equation}

\noindent where 

\begin{equation}
\begin{array}{c}
F = \vec{\sigma_{U}}(U, V) \cdot \vec{\sigma_{U}}(U, V)  \\
G = \vec{\sigma_{U}}(U, V) \cdot \vec{\sigma_{V}}(U, V)  \\
K = \vec{\sigma_{V}}(U, V) \cdot \vec{\sigma_{V}}(U, V)  \\
\end{array}
\label{metric_parametric}
\end{equation}

\noindent in which $\vec{\sigma_{U}}(U, V)$ and $\vec{\sigma_{V}}(U, V)$ are the gradient vectors at a point on the surface in parametric directions $U$ and $V$. From Equation \ref{EQ_h3d_to_h2d}, $E_{m_{AB}}$, $F_{m_{AB}}$, and $G_{m_{AB}}$ are obtained at point $m_{AB}$ of Figure \ref{measurements_in_2d}. Vector  $\vec{N}_{2D}$, normal to $AB$ (in 3D space), is computed using the metrics from Equation \ref{metric_parametric}: 

\begin{equation}
\vec{N}_{2D} = 
\left \{ 
\begin{array}{c}
U_{N_{2D}}  \\
V_{N_{2D}}  \\
\end{array}
\right \} 
 = \left[
 \begin{array}{lr}
G_{m_{AB}} & K_{m_{AB}}  \\
-F_{m_{AB}} & -G_{m_{AB}}  \\
\end{array}
\right]
 \cdot \vec{V}_{AB}
\label{normal_vector}
\end{equation}

\noindent in which $\vec{V}_{AB}$ is a unit vector from point $A$ to point $B$.


The input data for the present surface mesh generation algorithm is a parametric description of the surface and a polygonal description of the boundary of the surface patch to be meshed. The surface meshing algorithm has the following steps \cite{Miranda2002}: (a) the background quadtree, generated in Section~\ref{sec:SurfacesOctree}, is used to store local surface distortion metrics; and (b) a two-pass advancing front procedure is used to generate elements. In the first pass, the elements are generated based on geometric criteria, which produce well-shaped elements. In this phase, the background octree is also used to develop local guidelines for the new element sizes. Therefore, element sizes will take discretization error information into account. In the second pass, the elements are generated to guarantee a valid topological mesh, but with no good-shape requirements. For well-conditioned surface geometries, the second pass is usually not necessary.

A smoothing technique is used to improve the quality of the mesh by relocating nodes within a patch. A general formulation for this technique is based on a generic form of a weighted Laplacian function \cite{Foley1989}. 

The smoothing procedure is performed to satisfy the principle of equidistribution of nodes through the migration of nodes within a patch, which means that the nodes are relocated until they find an optimal position considering their neighborhood. Thus, relocating a particular node has an impact on its neighbors, and then the process should be repeated a number of times for all nodes until the optimal position for them is found. However, it is necessary to balance the quality improvement and the computational cost of the relocation. In this work, the smoothing procedure is repeated five times for all internal nodes, in order to achieve good quality results without impacting the performance. In addition, smoothing is done in the parametric space but takes into account the metrics of size distortion the between the parametric and 3D space.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Volume Refinement}\label{sec:VolumeRefinement}

3D mesh generation in each closed region of the model is based on a technique presented by Cavalcante-Neto \cite{Cavalcante-Neto2001} and is used to obtain tetrahedral elements in arbitrary domains. Similarly to the procedure applied to generate surface meshes, this one is based on an advancing front technique coupled to a recursive spatial decomposition technique (octree). Originally, the algorithm employed an independent background octree in each 3D region to control the distribution of the node points generated in the interior. In the adaptive methodology proposed here, the global background octree is used for this purpose.

The algorithm was designed to meet four specific requirements: to avoid producing elements with poor aspect ratios; to generate meshes conforming to existing triangular meshes at the boundary of a domain; to generate meshes exhibiting good transitions between regions of different element sizes; and to work properly for cases in which distinct boundary nodes are geometrically coincident (e.g. nodes on opposite faces of a crack).

The input to the algorithm is a triangular surface mesh, which describes the domain to be meshed. This mesh is obtained from the surface meshes on the boundary of a 3D region to be meshed. The algorithm steps are as follows:

\begin{itemize}
\item[$\bullet$] A two-pass advancing front procedure is applied to generate elements. In the first pass, elements are generated based on geometric criteria, producing well-shaped elements. The background octree presented in Section 4 is used to control the sizes of the elements and the position of the interior nodes. 
The octree determines an ideal position for an optimal node to form a new element. This ideal position defines a search region where an optimal node for the new element may be located. This region is a sector of a sphere whose center is the ideal position and whose radius is proportional to the octree cell size. If one or more existing nodes are inside this region, they are ranked based on a solid angle criterion, in order to get the best node for the new element. However, if no existing node is found, a new node is inserted at the ideal position and an element is generated using this node.  In the second pass, elements are generated based only on the criterion that they have valid topology. Here, any existing node that forms a valid new element can be used, regardless of whether it is close to the ideal position or not. However, the same quality criterion is used and the node that forms the best solid angle is chosen for the generation of the new element.
\item[$\bullet$] If the advancing front procedure cannot progress, a back-tracking strategy \cite{Cavalcante-Neto2005} is employed to delete some elements, and the procedure is restarted. It consists basically of back tracking a few steps in the mesh generation and deleting faces that hinder the front from converging. This creates better regions where valid elements can be then generated.

It is possible that the process of finding better regions may fail, for instance, if faces to be removed are part of the original boundary. When this occurs, other elements are deleted instead and the procedure is restarted. If a mesh still cannot be generated for this region, the algorithm fails and terminates. In principle, it is possible to create a boundary input mesh that forces the failure of the volume mesh generation. Such failure, however, has not yet been observed in ``non-contrived'' input, i.e., in any realistic input boundary meshes, such as the ones showed in the examples and many more tested so far.
\item[$\bullet$] Once a valid mesh is created, the quality of the element shapes is improved by using the standard Laplacian smoothing technique and locally deleting poorly shaped elements and those adjacent to them. The boundary contraction is then restarted.
\end{itemize}

After the generation of volumetric elements in all regions of the model as exposed above, a new error analysis is performed to assess the quality of the results. If convergence is not obtained, the whole adaptive process is repeated as described in the previous sections. The next section provides some examples of the proposed adaptive refinement process.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Examples}\label{sec:Examples}

This section presents examples of adaptive 3D finite element meshes that were generated using the proposed adaptive methodology. The main objectives of this section are to (1) analyze the visual aspect of the meshes generated, comparing the results with and without adaptivity; (2) compare the results of the convergence rates with results presented in the literature; and (3) comment on the obtained results. It is important to emphasize that this paper does not aim to compute the performance of mesh generators (surface and volume) or assess the quality of the elements generated, because these tasks were covered in previous works \cite{Miranda2002,Cavalcante-Neto2001}.

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-13.png}}
\caption{Example 1: Short cantilever under end shear.}
\label{FIG_exemplo_01}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2.5in]{./figure/Figure-14.png}}
\caption{Example 2: L-shaped domain under horizontal uniform force (Face y = 0).}
\label{FIG_exemplo_02}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2.5in]{./figure/Figure-15.png}}
\caption{Example 3: Biaxial bending of a column footing.}
\label{FIG_exemplo_03}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-16.png}}
\caption{Example 4: 3D Frame.}
\label{FIG_exemplo_04}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=1.8in]{./figure/rocker.png}}
\caption{Example 5: Bike suspension rocker.}
\label{FIG_exemplo_05}
\end{figure}

The adaptive strategy proposed in this paper results from the application of unstructured mesh generation techniques in surfaces and regions, combined with numerical errors associated with discretization. Numerical error estimators are implemented based on procedures developed for two-dimensional models \cite{Paulino1999} extended to three dimensions. These error estimators are supported by error estimation techniques widely adopted in the literature, called Superconvergent Path Recovery technique (SPR) \cite{Zienkiewicz1992a, Zienkiewicz1992b} and Recovery by Equilibrium in Patches (REP) \cite{Boroomand1997}. The numerical error estimators were implemented in a finite element numerical analysis program \cite{Martha2002}. It is worth stressing that these error estimation techniques can be easily and directly replaced by any other technique that is more recent or efficient, as this is supported by the object-oriented organization of the analysis program. In the following examples, SPR is employed for error estimation.

Figures \ref{FIG_exemplo_01} to \ref{FIG_exemplo_05} show five model examples to be tested employing the proposed adaptive strategy. The first three examples were obtained from the literature \cite{Le1997b}, and the adaptive process was applied by inserting additional nodes at the midpoint of the longest or quasi-longest line segment of the mesh, which bisects the original edges and generates new elements. In these examples, the adaptive refinement was carried out until the points were sufficient to allow comparison with reference results. The last two examples were chosen from real engineering problems, and the adaptive refinement was carried out until the relative error (3\%) was reached. The model examples are described below:

\begin{enumerate}
\item[1] \textit{Short cantilever under end shear}. A short cantilever beam is loaded in one extremity under uniform shear forces, while the other extremity is free. Due to symmetry, only half of the model is created, and the displacement in the $x$ direction along the face OHGA is set to zero; 
\item[2] \textit{L-shaped domain under tension}. An L-shaped domain is loaded under tension at the left extremity and free at the right extremity. Because of the distribution of internal stress, this model presents a singularity when it creates the L-form. Again, only half of the model is created due to symmetry, and the displacement in the face $x = 0$ is set to zero in the $x$ direction;
\item[3] \textit{Biaxial bending of a column footing}. This example is a column footing loaded at the top under uniform shear, face $z = 3$, and near the base under uniform pressure, face $x = 1.5$. All displacements of the footing face $z = 0$ are set to zero;
\item[4] \textit{Three-dimensional frame}. This example is composed by three regions with different materials, named Regions 1, 2, and 3. In all of them the cross section is 1 x 1. In Region 3, the model is loaded under uniform tension in the $x$ direction in the face $x = 6$. The displacements are fixed at one extremity of Region 1 in the face $y = 0$. This model covers a combination of external loading of tension (all regions), bending (Regions 1 and 2) and torsion (Region 1);
\item[5] \textit{Bike suspension rocker} \cite{JeanCuilliere}. This example is composed by four cylinders that are connected by a central body. Uniformly distributed forces, of unitary intensity in the $y$ direction are applied to the internal faces of the two top cylinders. Displacement restrictions $(u, v, w)$ are applied to the internal face of the lower cylinder.  
\end{enumerate}


\begin{table}
\centering
\caption{Mesh and convergence information for examples.}
\label{tab_mesh}       % Give a unique label
\begin{tabular}{p{1cm} p{1cm} p{1cm} p{1.0cm} p{1.0cm} p{1.0cm} }
\hline\hline\noalign{\smallskip}
Ex.--Step & \# of nodes & \# of elements & $\eta_{es}$ & $\xi_{avg}$ & $\xi_{dev}$  \\
\noalign{\smallskip}\hline\hline\noalign{\smallskip}
1--0 & 463 & 268 &  6.55 & 0.946 & 0.908 \\ 
1--1 & 4781 & 2785 &  2.09 & 0.197 & 0.369 \\
1--2 & 10479 & 6379 &  1.29 & 0.140 & 0.216 \\
1--3 & 19994 & 12376 &  0.85 & 0.105 & 0.134 \\
1--4 & 23835 & 14678 &  0.74 & 0.099 & 0.110 \\
\noalign{\smallskip}\hline\noalign{\smallskip}
2-0 & 268 & 137 &  7.97 & 1.170 & 1.387 \\
2-1 & 5760 & 3721 &  2.86 & 0.217 & 0.599 \\
2-2 & 12832 & 8525 &  1.45 & 0.155 & 0.292 \\
2-3 & 29251 & 19815 &  0.86 & 0.115 & 0.208 \\
\noalign{\smallskip}\hline\noalign{\smallskip}
3-0 & 1302 & 721 & 19.56 & 2.018 & 3.355 \\
3-1 & 16672 & 10682 &  6.82 & 0.710 & 1.164 \\
3-2 & 96537 & 65988 &  3.21 & 0.320 & 0.557 \\
\noalign{\smallskip}\hline\noalign{\smallskip}
4-0 & 229 & 112 & 13.71 & 2.102 & 1.768 \\
4-1 & 5419 & 3004 &  7.40 & 0.855 & 1.208 \\
4-2 & 57300 & 36996 & 1.35 & 0.177 & 0.849 \\
\noalign{\smallskip}\hline\noalign{\smallskip}
5-0 & 10303 & 6223 & 17.14 & 1.607 & 3.027 \\
5-1 & 273633 & 186238 &  3.97 & 0.307 & 0.733 \\
\hline\hline\noalign{\smallskip}
\end{tabular}
\end{table}


Figures \ref{FIG_adaptive_meshes_ex_01} to \ref{FIG_adaptive_meshes_ex_05} show the meshes generated by the adaptive refinement, considering quadratic tetrahedral elements. Similarly to the reference models Figures \ref{FIG_adaptive_meshes_ex_01}
 to \ref{FIG_adaptive_meshes_ex_03}, the meshes refine regions with singularities and regions with a high stress gradient. In Example 1 (Figure  \ref{FIG_adaptive_meshes_ex_01}), the mesh is refined along lines AB, BC, and CO, with singularities. In Example 2 (Figure  \ref{FIG_adaptive_meshes_ex_02}), the mesh is refined along the re-entrant corner line (AB), where the singularities is located, and coarsened on the face y = 1, where the constant tension stress is dominant. In Example 3 (Figure  \ref{FIG_adaptive_meshes_ex_03}), the mesh is refined mostly at the junction between the column and the footing (lines CD, DF, EF, and CE), where singularities are located. Other refinement locations are the top of the column under uniform shear and near the base under uniform pressure, creating a singularity along line AB. In Example 4 (Figure \ref{FIG_adaptive_meshes_ex_04}), the adaptive refinement presents an important characteristic of the proposed methodology: to coarsen locations of low gradient stress. In Region 3, the mesh is refined in the first step and coarsened in the second (final) step, because the constant tension stress is dominant. In this example, only one background octree was used to hold the new sizes of the elements, another feature of the proposed methodology. The final mesh is coherent with the FE analysis, because the mesh is more refined in Region 1, where complex stress fields are present. In example 5 (Figure \ref{FIG_adaptive_meshes_ex_05}), the mesh is refined in the intermediate cylinder and only part of the lower cylinder, where there are concentrations of stress. This example also demonstrates the importance of considering the curvatures of the supporting surface in the adaptive refinement, which is another characteristic of the present methodology. Figure \ref{FIG_adaptive_meshes_ex_05_ref} shows a detail of the mesh refinement, where part (on the bottom) of the lower cylinder is refined based on curvature.

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-18.png}}
\caption{Adaptive refinement meshes for Example 1}
\label{FIG_adaptive_meshes_ex_01}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-19.png}}
\caption{Adaptive refinement meshes for Example 2}
\label{FIG_adaptive_meshes_ex_02}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-20.png}}
\caption{Adaptive refinement meshes for Example 3}
\label{FIG_adaptive_meshes_ex_03}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/Figure-21.png}}
\caption{Adaptive refinement meshes for Example 4}
\label{FIG_adaptive_meshes_ex_04}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/rocker_exem.png}}
\caption{Adaptive refinement meshes for Example 5}
\label{FIG_adaptive_meshes_ex_05}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=2.0in]{./figure/detalhe_rocker.png}}
\caption{Detail of mesh refinement for Example 5.}
\label{FIG_adaptive_meshes_ex_05_ref}
\end{figure}


Table \ref{tab_mesh} shows mesh and convergence information for all examples. ``Ex.-Step'' means example and current step of adaptivity, $\eta_es$ is the estimated relative percentage error, $\xi_{avg}$ and  $\xi_{dev}$ are estimated and standard deviation error ratio of finite element model, respectively. Note that all convergence information decreases to a minimum error. 

Figure \ref{FIG_rate} shows the results of convergence rates compared to those of the reference \cite{Le1999}. All results obtained herein converge faster to a small relative error. For instance, in Example 1, the reference paper required three steps to achieve a relative error close to 1\%, while the present work needed only two steps with adaptive refinement. The difference between these results is due to the fact that the present work employs a new global generation of tetrahedral elements in all regions, allowing the new elements to have greater of shape and size flexibility, whereas using a local refining approach, as the reference work did, the shape and size of the new elements depend on the old mesh. All examples presented herein achieved a relative error close to 2$\sim$3\% with only two adaptive refinements.

\begin{figure}[h]
\centering
\subfigure[FIG_rate_ex_01][Example 1]{\includegraphics[width=3.30in]{./figure/Figure-24.png}}
\subfigure[FIG_rate_ex_02][Example 2]{\includegraphics[width=3.30in]{./figure/Figure-25.png}}
\subfigure[FIG_rate_ex_03][Example 3]{\includegraphics[width=3.30in]{./figure/Figure-26.png}}
\caption{Comparison of convergence rate.}
\label{FIG_rate}
\end{figure}

\begin{figure} 
\centerline{\includegraphics[width=3.34in]{./figure/JMesh_time.png}}
\caption{Time performance of 3D mesh generator.}
\label{FIG_time_Jmesh}
\end{figure}

To estimate the performance of the proposed methodology, Figure \ref{FIG_time_Jmesh} shows the elapsed processing times of 3D mesh generator for the examples. Here only the volumetric refinement step is included, because this is the step where it takes longer to generate the final finite element model. These processing times were measured on a Pentium Core i7-1.73GHz PC with 6 GB of RAM, under Windows 7 operating system. Note that the graph is in $log-log$ scale. Using this graph, it is possible to estimate the time required for 1 million elements, for example, which would take nearly 17 minutes. However, it is better to generate 10 times 100 thousand elements in approximately 7 minutes.       


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}\label{sec:Conclusion}

This paper described a methodology for adaptive generation of three-dimensional finite element meshes, using geometric modeling with multi-regions and parametric surfaces. The mesh adaptive process involves three steps: (1) FE analysis with error estimation; (2) construction of a structure to store the new sizes of the FE; and (3) refinement of the FE model. The approach adopted is the global refinement of the whole model in each adaptive refinement, using a background octree structure. After the construction of the octree, the new model is geometrically re-discretized employing a hierarchical curve, surface and volume refinement.

Some examples have demonstrated important characteristics of the proposed methodology, namely:

\begin{itemize}
\item[$\bullet$] The ability to refine and coarsen in regions of high and low response gradients;
\item[$\bullet$] The use of only one background octree for all regions of the model, allowing a smooth transition between regions and elements;
\item[$\bullet$] The hierarchical refinement of curves, surfaces, and volumes;
\item[$\bullet$] The consideration of curvatures of the supporting curves and surfaces in the adaptive refinement.
\end{itemize}

In addition, comparing the results obtained by the present work with the ones in the referred literature, the former converge faster to a lower relative error, because the 3D mesh generator used in this work has more freedom to create new elements based on desired element sizes. This characteristic generates a desired mesh with the application of only a few steps of the adaptive refinement. In contrast, using a local refinement strategy, (local) element manipulations restrict the shape quality of new elements.

Obviously, the current 3D mesh generation takes more time to create new elements, because the whole FE model must be created at each step. In this work, most of meshes were generated in less than one minute of clock time. However, in models that require a large number of elements (one million or more), time consumption increases exponentially, as shown in Figure \ref{FIG_time_Jmesh}. This is a significant limitation of the global approach that can be solved in two manners: (1) decomposing the domains into sub-domains and applying the mesh generator to each sub-domain; or (2) using a parallel 3D mesh generator. The latter option is our future work and is currently under development.          



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgment}
The authors would like to thank the National Council for Scientific and Technological Development (CNPq), the CAPES Foundation of the Ministry of Education of Brazil, the Computer Graphics Technology Group (Tecgraf), and the Pontifical Catholic University of Rio de Janeiro (PUC-Rio) for the financial support and for providing the necessary space and resources used during the development of this work. We thank ANP, Petrobras, and the Laboratory of Scientific Computing and Visualization (LCCV) of the Federal University of Alagoas for granting access to the GradeBR/UFAL cluster computing resources of Rede Galileu.
\end{acknowledgment}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here's where you specify the bibliography style file.
% The full file name for the bibliography style file 
% used for an ASME paper is asmems4.bst.

\bibliographystyle{asmems4}
\bibliography{miranda_etal}

\end{document}
